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Abstract 

O The Polyakov loop extended Nambu and Jona-Lasinio model (PNJL model) in a 

mean field framework shows astonishingly good agreement with lattice QCD calcu- 
qj • lations which needs to be better understood. The present work reports on further 

developments concerning both Polyakov loop and mesonic fluctuations beyond mean 
field approximation. Corrections beyond mean field are of special interest for the 
thermal expectation values of the Polyakov loop ($) and its conjugate ($*), which 
^ • differ once the quark chemical potential is non-zero. Mesonic fluctuations are also 

\ considered with emphasis on the role of pionic modes. 

in 

1 Introduction 

(N 

Exploring the thermodynamic properties of strongly interacting matter has become a 
central theme of high-energy nuclear physics in recent years. On the theoretical side, 
great progress has been achieved thanks to lattice calculations solving discretised QCD 
numerically. At the present stage a major part of the numerical expense at finite quark 
chemical potential is caused by the fermion sign problem. The three most promising 
ways to address this difficulty are multi-parameter re-weighting techniques [1, 2], analytic 
continuation from imaginary chemical potentials [3, 4] and Taylor series expansion methods 
[5,6,8,9,11-13]. 

It is an important task of effective field theories and models to reveal principal mech- 
anisms and their functioning behind the otherwise hidden mechanisms of lattice QCD. 
Alongside with other approaches [14-16] the Nambu and Jona-Lasinio (NJL) model [17— 
24] is an approach that successfully describes spontaneous chiral symmetry breaking. We 
use the NJL model with Nf = 2 quark flavours as one of our starting points. In the NJL 
model gluonic degrees of freedom are "integrated out" . The role of the gluons is assumed 
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to be modelled in part by a local effective quark colour current interaction. From this ef- 
fective interaction a Fierz transformation generates various quark-antiquark and diquark 
coupling terms. By integrating out the gluons the local SU(3) C gauge symmetry is lost. 
As a consequence the NJL model, equipped only with a global SU(3) C symmetry, does not 
feature confinement. 

To bring aspects of confinement back into the model an additional homogeneous tem- 
poral background field with standard SU(3) C gauge invariant coupling to the quarks is 
introduced [17,18], implementing the Polyakov loop. This generalised NJL model, with 
Polyakov loop dynamics incorporated, is called PNJL model. In the limit of static quarks 
(i. e. in pure gauge QCD) the Polyakov loop serves as an order parameter for confinement. 
In this limit the Z(3) centre-symmetry of the SU(3) C gauge group is unbroken, and the 
deconfinement transition is connected with the spontaneous breakdown of this symmetry. 
In the presence of dynamical quarks the Z(3) centre-symmetry is broken explicitly, such 
that the deconfinement transition is no longer a phase transition in the strict sense. Never- 
theless, the Polyakov loop shows a rapid crossover near the deconfinement transition, still 
permitting to use the Polyakov loop as a measure for deconfinement. The confining gluon 
dynamics that was lost in the NJL model is now re-introduced via an effective potential. 
This potential is part of a Ginzburg-Landau model for confinement in the static quark 
limit. The information necessary in order to specify the effective potential is extracted 
from pure glue lattice QCD calculations [13]. 

In Sec. 2 we review the PNJL model [19, 20, 22-24] at the level of the mean field 
approximation. When implementing the Polyakov loop extension to the NJL model, it 
is important to pay special attention to the fermion sign problem. As the Polyakov loop 
is coupled to the NJL model in analogy to QCD via minimal substitution the fermion 
sign problem in QCD and in the PNJL model appears on equal footing. Issues arising 
from the fermion sign problem in the PNJL model are discussed in Sec. 3. In the present 
paper these issues are addressed more explicitly than in previous work [19,25,26]. The 
method developed in this improved treatment introduces a systematic expansion around 
a leading order (mean field) approximation defined such that all physical quantities are 
real in this limit. It is demonstrated that the dynamics beyond mean field can be treated 
perturbatively. A detailed derivation of these perturbative corrections is presented in the 
appendix. The perturbative method is applied in Sec. 4 to the PNJL model to investigate 
the effects of the complex phase of the action on different quantities. Here the expectation 
values of the Polyakov loop and its conjugate are of special interest. In the present 
analysis, the split of the expectation values of the Polyakov loop and its conjugate at 
non-zero chemical potential arises once fluctuations of the fields are taken into account. 
Quantities like susceptibilities, in which mean field contributions partially or completely 
cancel, are sensitive to the corrections beyond mean field, as we shall discuss. Finally 
in Sec. 5 further corrections beyond Hartree approximation are estimated, generated by 
propagating mesonic (quark-antiquark) modes. The lightest meson mode, the pseudoscalar 
pion with its approximate Nambu-Goldstone character, is the leading correction in this 
sector. Sec. 6 presents our conclusions and an outlook. 
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2 The PNJL model 



The two-flavour PNJL model including diquark degrees of freedom [25] is derived from 
the Euclidean action 

S E (4>,^,<t>) = J dr J d 3 x ^d T ^ + H(ip,^,<f>)\ +5S E (<f>,T) (2.1) 

with the fermionic Hamiltonian density 1 : 

H = -# f (a • V + 74 m - A 4 ) ip + V(tp, ^) , (2.2) 

where ip is the Nf = 2 doublet quark field and tuq = diag(m M , m^) is the quark mass matrix. 
The quarks move in a background colour gauge field A4 = iAq, where Aq = 6^0 gAa t a with 
the SU(3) C gauge fields Aa and the generators t a = X a /2. The matrix valued, constant 
field A4 relates to the (traced) Polyakov loop as follows: 

$ = -ptr c L with L = exp (i [ drA^) and j3 = ^ . (2.3) 
A c \ Jo J T 

In a convenient gauge (the so-called Polyakov gauge), the matrix L is given a diagonal 
representation 

L = exp [i (03 A 3 + (f>8 A 8 )] • (2.4) 

The dimensionless effective fields ^3 and (ps introduced here are identified with the Eu- 
clidean gauge fields in temporal direction, divided by temperature: $3 = A^ /T and 

(8) 

4>s = A 4 /T. These two fields are a parametrisation of the diagonal elements of SU(3) C . 
As such the "angles" ^3 and </>s necessarily have to be real quantities in order to sustain 
the unitarity of the group. An alternative parametrisation of the diagonal elements of 
SU(3) C is given by the Polyakov loop, <E> = |tr c L, and its conjugate, <I>* = |tr c LL 

The piece 5Se = of the action (2.1) carries information about the gluon dynam- 
ics. The potential 1A effectively models the confinement-deconfinement transition and the 
region up to temperatures of roughly T < 2 T c in quarkless, pure gauge QCD on the mean 
field Ginzburg-Landau level. At temperatures very far above the transition a description 
of the thermodynamics with just the two order parameters $ and is not appropriate 
as transverse gluons will become important. Transverse gluon degrees of freedom cannot 
be described by Polyakov loops. 

The Polyakov loop is an order parameter for confinement in SU(3) gauge theory. In the 
confined low temperature phase the expectation value of the Polyakov loop vanishes, ($) = 
0, while (<J>) 7^ implies deconfinement. Let To be the critical temperature separating the 
two phases. As previously mentioned the symmetry which is restored at T < To and 
broken above To is the Z(3) centre-symmetry of SU(3) 2 . 

Therefore, the Landau effective potential describing the dynamics, the Polyakov loop 
potential U(&,T), has to be Z(3)-symmetric in <I>. The basic building blocks for such a 
potential are $*<I ) , $ 3 and < I>* 3 terms. The potential used here differs from the simplest 



1 a = 70 7 and 74 = 470 in terms of the standard Dirac 7 matrices. 

2 The centre of SU(3) contains all those SU(3) elements that commute with all other SU(3) elements, 
i. e. the elements e l ~ h 1, with fceZ constituting a Z(3) subgroup of SU(3). 
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ansatz generating a first order phase transition as it is implemented in [19]. Instead we 
use the ansatz given in [25, 26] motivated by the SU(3) Haar measure: 

^l^lll = -l a (T)^*^ + b(T) ln[l-6$*$ + 4(V 3 + $ 3 ) -3($*$) 2 ] , (2.5) 
where the temperature dependent prefactors are given by 

a(T) =a + ai (^j + a 2 (^j ' and b(T) = b 3 (^j ° . (2.6) 

The logarithmic divergence near $ — > 1 properly constrains the Polyakov loop to 
values attainable by the normalised trace of an element of SU(3). The parameters of 
U(§, T) are chosen such that the critical temperature of the first order transition is 
indeed equal to T (fixed at 270 MeV [29]) and that , $ -> 1 as T -> oo. 

The numerical values using these constraints are taken as given in Refs. [25, 26] 

a = 3.51 , ai = -2.47 , a 2 = 15.2 , b 3 = -1.75 . 

The resulting uncertainties are estimated to be about 6% for ai, less than 3% for a 2 
and 2% for 63. The value ao = chosen here reproduces the Stefan-Boltzmann limit. 
This is not mandatory, of course, since the high-temperature limit is governed by (trans- 
verse) gluonic degrees of freedom not covered by the Polyakov loop which represents the 
longitudinal gauge field. Alternative parametrisations of U are possible, such as the two- 
parameter form guided by the strong-coupling approach [18], which has a different high 
temperature limit. In the present context these differences are not crucial as we system- 
atically restrict ourselves to temperatures close to the transition region, T < 2T C , where 
different forms of IA give remarkably similar results as pointed out in Ref. [27]. 

In Fig. 1 we plot the Polyakov loop potential using the parametrisation given in 
Refs. [25,26] at T = T = 0.27 GeV. This illustrates the Z(3) symmetry. The single 
minimum at T < Tq becomes degenerate with three minima at T = To. Above To only 
these three minima survive. Of course upon spontaneous breakdown of the Z(3) centre- 
symmetry, the three minima and the Z(3) centre-symmetry of the potential remain intact 
even though the vacuum expectation value does no longer show the symmetry of the 
potential. 

The NJL interaction term V in Eq. (2.2) includes chiral SU(2) x SU(2) invariant 
four-point couplings of the quarks acting in pseudoscalar-isovector /scalar- isoscalar quark- 
antiquark and scalar diquark channels: 



(V^) 2 + (^75rV) 2 ] - y [(^C 75 r 2 A 2 ^ T ) ^ T 75r 2 A 2 C^] , (2.7) 



where C is the charge conjugation operator. These interaction terms in Eq. (2.7) are 
obtained from a local colour current-current interaction between quarks, 

Ant = -G c (V^ a V0(V ; 7^>) , 

by a Fierz transformation which relates the coupling strengths G and H as G = which 
we choose not to alter 3 . 



3 Additional terms generated by the Fierz transformation are of no importance in the present context 
and will be omitted. 
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The NJL model with two quark flavours is usually modelled with three parameters, a 
current quark mass m U; d, a local four quark coupling strength G and a three-momentum 
cutoff A. The parameters used here are the ones used in [19, 25, 26]: 

m Ujd = 5.5MeV, G = ^-H = 10.1 GeV~ 2 , A = 0.65GeV, 

fixed to reproduce the pion mass and decay constant in vacuum and the chiral condensate 
as m„ = 139.3 MeV, f w = 92.3 MeV and (^ u tp u ) = -(251 MeV) 3 . 

To evaluate the thermodynamic properties of the model the quark degrees of free- 
dom are integrated out. New auxiliary fields are introduced by bosonisation, absorbing 
quark-antiquark and quark-quark (antiquark-antiquark) correlations. These are a scalar- 
pseudoscalar field (a, tt) and a diquark (antidiquark) field A (A*). The resulting thermo- 
dynamic potential then reads 

n = ^S hos =U(<£, $*,T)-^ J ^Trln[pS- 1 (^n,p)}+^ + ^ , (2.8) 

where the Matsubara sum runs over uj n = (2n + 1)ttT reproducing antiperiodic boundary 
conditions in the Euclidean time direction. The inverse Nambu-Gor'kov propagator S^ 1 
in Eq. (2.8) is defined by 



S 1 (iuj n ,p) 



«7o w„ - 7 -p - m + 70 {n - iA 4 ) A75T2A2 

-A*7 5 t 2 A 2 no w n - 7 • P ~ m ~ 7o (m - iM) 

(2-9) 

The mass of the quark-quasiparticles is given as in the standard NJL model by the gap 
equation 

m = m - ( a ) = m - G^ip) . (2.10) 

The Matsubara sum is evaluated analytically. The quasiparticle energies emerging in 
this procedure are related to the solutions of det [5 _1 (po)] = 0. The bosonised action 
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then reads 



T a 2 A* A 

n = v Sbos = U T) + — + — 



with six distinct quasiparticle energies 

Ei,2 = e{p) ± h , 



E 3A = V(e(p)+Ar) 2 + |A| 2 ±iT0 3 , 

£5,6 = V(e(p)-Ar) 2 + |A| 2 ±iT0 3 , (2.12) 
where e{p) = \/ p 2 + m 2 . Additionally we have introduced 

£ 6 = /i + 2iri|, /ir = /i-iT^. (2.13) 

The energy difference A£j- is defined as the difference of the quasiparticle energy and the 

energy of a free fermion, £0 = \J p 2 + m 2 ,: AEj = Ej — Eo ± /i. The form of the bosonised 
action, Eq. (2.11), does not allow to factor out the Polyakov loop fields $ and $*, as it was 
done in Ref. [19]. Instead we keep the form of Eq. (2.11) using ^3 and 4>g with ^3, </>g 6 R. 4 

The introduction of the Polyakov loop outlined above formally leads to a complex 
valued action as soon as fi 7^ 0. This phenomenon is usually called fermion sign problem. 
Due to the connection of ^3 and (ps to the QCD colour gauge group SU(3) C , we must 
require (ps and <fis to be real fields at all times. The mean field approximation, ^mf, of 
the thermodynamic potential f2 must be introduced such that it satisfies this constraint 
imposed by the gauge group. Identification up to a constant of Omf with the (real) 
pressure p in this approximation then requires that non-Hermitian structures of the inverse 
quasiparticle quark propagator do not contribute. One way to establish such a lowest order 
approximation is to use the real part of the thermodynamic potential in the mean field 
equations. 

The necessary condition for the minimisation of the effective action in a standard 
situation is, in general, 

where d{ stands for the fields representing the relevant degrees of freedom (in our case: 
9 = (a, A, 03, (f)s))- In order to always comply with ^3, 0g G R we define the mean field 
thermodynamic potential, with Qq of Eq. (2.11), by 



f^MF = Re[Oo] = Re 



y <5t>OS 



(2.15) 



As the parameter space of 4>z and (/>$ is periodic there are different parameter sets representing the same 
physics. We use the (triangular shaped) domain {(<j> 8 ^ —J.)A(<f> 8 ^ VS^a + ^^A^s ^ \/3(-<?(>3 + 4f ))}. 

Note that the periodic domain of L and iJ is 3!-times larger than the domains for $ and $* (or equivalently 
4>3 and 4>g) due to the trace's invariance under unitary transformations of L. 
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0.0 0.5 1.0 1.5 

T/T c 

Figure 2: Chiral condensate normalised to its value at temperature T = (dash-double- 
dotted line) in the NJL model with massless quarks, and Polyakov loop (<3?) in the pure 
gauge model (dashed line). The PNJL model (with non-zero quark masses) shows dynam- 
ical entanglement of the chiral (solid line) and Polyakov loop (dash-dotted line) crossover 
transitions. For comparison lattice data for the Polyakov loop in pure gauge and full QCD 
(including quarks) are also shown [30] . 



The mean field equations then read 

d^MF d Re [Op] 

d (a, A, 03, 0s) d (a, A, </> 3 , 4> 8 ) ' 1 " ' 

The hereby neglected imaginary part of this derivative, gl^~£^~h ; W1U be taken into 
account by writing down a series in powers of this residual gradient. In addition it is also 
possible to correct for deviations of the potential from a gaussian shape which is assumed 
for the mean field approximation. As explained in the appendix it is most convenient to 
consider both types of corrections simultaneously using Feynman graphs to construct all 
possible terms. 

A variety of PNJL model results (equations of state, phase diagrams, susceptibilities) 
have been obtained in previous calculations [25, 26, 31] based on the mean field equations 
(2.16). At this point it is instructive to examine how chiral and Polyakov loop dynamics 
cooperate to produce crossover transitions (at zero chemical potential) which end up in 
a narrow overlapping range of temperatures (see Fig. 2). In isolation, the pure gauge 
Polyakov loop sector and the NJL sector in the chiral limit show first (second) order phase 
transitions with critical temperatures far separated, as demonstrated by the dashed (dash- 
double dotted) lines in Fig. 2. When entangled in the PNJL model, these transitions (with 
non-zero quark masses) move together to form a joint crossover pattern. 

3 Fluctuations and corrections beyond mean field 

The discussion of fluctuation corrections to the mean field approximation has two parts: 
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• corrections arising from the imaginary part of f^o an d involving the Polyakov loop 
<I> and its complex conjugate <£*; 

• corrections from dynamical fluctuations involving propagating meson fields, with 
emphasis on the pion. 

The first item is the primary topic of this present section. The second item will be relegated 
to a separate Section 5. 

The mean field equations (2.16) establish a leading order approximation satisfying the 
reality constraints on the fields </>3 and <ps- Eq. (2.16) is at the same time the necessary 
condition for the maximisation of the modulus |e E | of the thermodynamic weight in the 
path integral 5 : 

The direct connection $ = <J?(</>3, 0g) of the two parametrisations, (fa, fa) on the 
one hand and ($, $*) on the other hand, is lost once we step away from mean field and 
calculate thermodynamic expectation values: ($) 7^ Q((fa) , (0s)). 6 This observation 
is crucial when comparing the present method of approximation to schemes in previous 
publications [19,20,22-24]. In these publications the fields fa, fa have been replaced by 
<£, before doing mean field approximation. This implies that (<!>) and (<£*) (and not 
<£, <£*) are treated as independent mean field degrees of freedom. The minimisation of f^o 
is then performed requiring that (<3?) and (<£*) are real quantities. In such approximation 
schemes it is therefore not possible to find a way back to the (real) quantities ^3, fa: ($) 
and ($*) already comprise fluctuations of fa, fa £ IR. In other words, the definition of 
the lowest order approximation (which is usually referred to as mean field approximation) 
is different in Refs. [19, 20, 22-24] and this work. The definition of the lowest order (mean 
field) approximation in this work allows to strictly separate contributions originating in 
constant and fluctuating parts of the fields. 

The combination of the constraints G H and fa, fa £ R allow only certain limited 
configurations of fa and fa. In mean field approximation the condition fa, fa £ H implies 
that, due to Eq. (2.4), <I> and are the complex conjugates and we find ($} MF = $mf, 
(<3?*) MF = ^mf- At /k = the Polyakov loop <E> and its complex conjugate $* are treated 
equally due to charge conjugation invariance. It follows that <3?mf = $mf £ R in mean 
field approximation, fixing fa = 0. The Polyakov loop effective potential U = U(T, <3?, $*) 
in its parametrisation (2.5) is minimal for <I> = <E>* at fixed |$|. We find that the Polyakov 
loop potential U is always strong enough to keep <E> = <I>* or, equivalently, fa = 0. Not all 
parametrisations of IA will maintain this solution. If the curvature of the potential IA is 
not strong enough, the solution fa = becomes instable, and fa = is the position of a 
local maximum of the potential. In Ref. [20] a symptom of this fact has been described: 
the susceptibility associated with Re$ may become negative. The potential used in the 
present work does not show such deficiencies. 

After these preparatory remarks we proceed to develop a calculational scheme which 
systematically treats corrections to the mean field approximation, Eqs. (2.15) and (2.16). 
The basic idea is, as usual, to expand the thermodynamic potential around its mean 

5 In analogy to the procedure in Minkowskian space-time one might argue that the complex phase 

needs to become stationary. Taking the thermodynamic limit one observes that the stationary phase field 

configuration is favoured over any other configuration by the factor , while the absolute value is favoured 

v_ 

over other configurations by a factor e t . 

6 Recall that the thermal expectation value (• • •} is a weighted sum (integral) over various thermal field 
configurations. The equality (<!>} = $(((^3} , (</>s)) holds only if $ and $* are linear functions of (f>3 and 4>g. 
This is not the case. 
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field limit Omf- Technical details of the derivation are summarized in the appendix. The 
result including next-to-leading order is 



(dn \ T 


'd 2 n ~ 


1 dn 


{de J 


de 2 _ 


de 



(3.1) 



starting from the (complex) Qq of Eq. (2.11), with f^MF defined by Eq. (2.15), and with 
#mf determined by Eq. (2.16). The gradients dQo/de are understood with the set of 
field variables e = (9i) = (cr, A, 03, (/)%) arranged in vector form, and d 2 £lo/de 2 stands 
for the matrix (9 fio/9^9Sj). The correction term in (3.1) is taken using the mean field 
configuration, 6* = #mf- Note that this term takes care of the contributions from Imf2o in 
such a way that Q remains a real quantity. 

The thermal expectation value (/) of a physical quantity / is calculated according to 



(/) = /(0mf) - 



(dn \ T 


'd 2 n ~ 


1 Of 


{de J 


de 2 _ 


' de 



(3.2) 



Applications will now be given, in particular, for the Polyakov loop and for susceptibilities. 



4 Results 

The numerical calculations presented in this section are performed in the thermodynamic 
limit, i. e. in leading order of the S-expansion, and up to first order in the (^-expansion, 
as explained in detail in the appendix. The thermodynamic potential is determined by 
Eq. (3.1). Thermal expectation values are computed using Eq. (3.2). 

4.1 The Polyakov loop (<&) and its conjugate (<&*) 

With the mean field definition (2.15) the Polyakov loop expectation values (<5) and (<J>*) 
turn out to be equal in this limit, given the reality constraint on Omf- It is the corrections 
from Im f2o induced by the temporal gauge fields which cause the splitting of (<fr) and ($*). 

The difference ($>*} — ($) vanishes at zero quark chemical potential \i and has the same 
sign as fi, in agreement with results of Ref. [28]. As can be seen from Fig. 3 the difference 
(<J>*) — ($) is pronounced around the phase transitions. In the upper left panel of Fig. 3 
the influence of the first order phase transition separating the chiral and the diquark phase 
at low temperature can be seen as a jump in both (<J>) and (<£*). The second order phase 
transition separating the diquark regime from the high temperature quark-gluon phase 
can be identified as a kink in the lower right panel of Fig. 3. 

4.2 Susceptibilities and phase diagram 

A susceptibility Xg involving a quantity g is defined by 

X 2 g = V ((a - (g)) 2 ) = V (V) - (g) 2 ) . (4.1) 

Susceptibilities of special interest in the present context are the ones related to the dy- 
namical quark mass, m = tuq — a, and to the Polyakov loop. They are expressed in terms 
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t/t c t/t c 



Figure 3: Examples of thermal expectation values of the Polyakov loop ($} and its conju- 
gate (<£*). In the upper row ($) and ($*) are plotted as functions of the chemical potential 
H at constant temperature T. Below (<&) and ($*) are plotted as functions of temperature 
T at constant chemical potential fi. 
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of the inverse matrix of the second derivatives of the full thermodynamic potential f2: 



X 2 m = V (V) - (m) 2 ) = T dn 

x! = ^(<<& 2 >-<*> 2 ) = T 



dQidQi 

1 J J m,m 



d 2 n 



T 


r d 2 n ' 


-i r 


"<9 2 O fu if 


-l 


T 

+ 4 


r a 2 ^ " 


-l 


1 




,86*86 j_ 









(4.2) 
(4.3) 
(4.4) 



These susceptibilities are calculated using the graph rules given in Tab. 2 of the appendix 
which lead to the following explicit form: 



Xg 



= T 





'd 2 n ~ 


1 dg 




de 2 _ 


' de 



- 2T 






'd 2 n ~ 


1 dn \ 


( 


oe 2 _ 


de 

/ k 



(4.5) 



9=8 M f 



Here g stands for m or 3>, respectively. The first term in Eq. (4.5) is the susceptibility of 
the gaussian theory whereas the other terms are interpreted as corrections. 

The susceptibilities xm an d x$ serve as indicators for boundaries between phases 
when drawing a phase diagram in the plane of temperature and chemical potential. For 
smooth crossover transitions, such boundaries are not rigorously defined. Several criteria 
can be used to determine a transition line separating the region of spontaneously broken 
chiral symmetry from the quark-gluon phase. We use here the maxima of the chiral 
susceptibilities \M and of the Polyakov loop susceptibility XRc<i> in comparison with the 
maximal slopes dm/dT and d (<&* + <J>) /dT of the corresponding quantities which act, 
respectively, as order parameters in the limiting situations of exact chiral SU(2) x SU(2) 
symmetry or Z(3) symmetry. 

Fig. 4 shows a comparison of crossover transition lines found with the two criteria 
just mentioned. As both criteria are linked via the quadratic term in the action, all 
curves finally converge to the same point, the critical point (here in the absence of diquark 
condensation). A singularity in the second derivative of the action (or equivalently in 
the propagator) enforces this unique intersection point where the specific heat and other 
quantities show singular behaviour. 

Comparing our Fig. 4 with corresponding results in other publications (see Fig. 16 in 
Ref. [20] and Fig. 4 in Ref. [21]) one finds that the detailed behaviour of the deconfinement 
crossover transition line depends sensitively on the parameter choice and regularisation 
prescription. In the present case of a strong coupling a joint course of chiral and de- 
confinement crossover line is observed. When the coupling becomes weaker (e. g. due to 
parameter choice and regularisation prescription as in Refs. [20,21]) the transition lines 
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Figure 4: Comparison of the transition lines obtained by determination of the maximum 
in the chiral susceptibility (left panel solid line) and the Polyakov loop susceptibility XRe $ 
(right panel solid line) with the transition lines fixed by the maximal change with respect 
to temperature of constituent quark mass (left panel dashed line) and average of the real 
part of the Polyakov loop \ ($ + <&*) (right panel dashed line). 



may deviate. In any case one should note that such deviations appear at large chemical 
potentials approaching the typical cutoff scale of the model. Any conclusions drawn at 
such energy or momentum scales should be handled with care. 

In Fig. 5 we show the chiral and the Polyakov loop susceptibilities as functions of 
temperature at vanishing quark chemical potential (left panel) and compare them to the 
temperature derivatives of the constituent quark mass and the Polyakov loop (right panel). 
If we consider the behaviour of xm, XRe$ and x$ at T — > we find that XRc$ is finite, while 
Xm and x$ vanish. This can be explained by the fact that (Re <E>) 2 = \(^ 2 + 2 |<3?| 2 + $* 2 ) 
contains a U(l)-symmetric term |$|. As the U(l) symmetry incorporates the Z{3) centre 
of SU(3) C this term does not have to vanish once the Z{3) symmetry is fully restored at 
T = 0. 7 The width of the peak in the temperature derivative of the dynamical quark mass 
m = rriQ — a suggests that this crossover is influenced by the crossover of the Polyakov loop. 
At finite current quark mass mo the PNJL model produces an approximate coincidence of 
the peaks in the susceptibilities of the Polyakov loop and the constituent quark mass m, 
consistent with the pattern observed in Fig. 2. 

A comparison of the phase diagram obtained in mean field approximation in Ref. [25] 
and the phase diagram including corrections to the order (3 < 1 shown in Fig. 6, explicitly 
approves that corrections to the phase diagram due to the fermion sign problem are indeed 
small [25]: the influence of Imfio and the splitting of (<&*) and (<£) are rather modest. 

4.3 Moments of the pressure 

One benchmark for the PNJL model is its surprising capability of reproducing the trends 
of lattice QCD calculations. 8 One way to handle the fermion sign problem in lattice QCD 
is to expand the calculated pressure about \i = in a Taylor series. Such an expansion is 

7 The authors thank Chihiro Sasaki for pointing this out to them. 

s Note however the discussion concerning the dependence on quark masses in Ref. [31]. 
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Figure 5: The chiral susceptibility xm (left panel solid line) and the Polyakov loop suscep- 
tibilities XRe$ (left panel dashed line) and x$ (left panel dotted line) plotted as functions 
of temperature at vanishing quark chemical potential. These susceptibilities defined by 
Eqs. (4.2, 4.3, 4.4) and evaluated using Eq. (4.5) are compared here to the derivative of 
the constituent quark mass (right panel solid line) and the expectation value of the real 
part of the Polyakov loop (right panel dashed line) with respect to temperature. 
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Figure 6: Phase diagram implementing corrections to the order (3 < 1. Solid lines: 
crossover transition of the susceptibility related to the real part of the Polyakov loop, 
dashed lines: first order phase transition, and dotted: second order phase transitions. 
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Figure 7: The moments of the pressure with respect to ^ as defined in Eq. (4.6). C2 is 
shown in the left panel, C4 is displayed to the right where cq is shown in the inset. The 
data deduced from lattice computations are taken from [9]. 



given in Ref. [9]: 



- ^ c n (T) {-) with c n (T) - - — — - 

with even n as the situation is charge conjugation invariant. Specifically: 



(4.6) 



1 d 2 {p/T A ) 
c 2 = 77 



2 0(/VT)2 

1 a 6 (p/r 



1 a 4 (p/r 4 ) 

c 4 = 7T7 



/u=0 



720 9(/x/T) 6 



24 d(fi/T) 4 
1 9 8 (p/T 



At=0 



/i=0 



40320 a(^/T) 8 



(4.7) 



H=0 



The pressure in the PNJL model is evaluated by subtracting the divergent vacuum con- 
tributions of the thermodynamic potential: 



p = -(n-n(r = o)) 



(4.8) 



Results for C2, C4 and C6 are shown in Fig. 7. In comparison with the plots for c n 
presented in a previous paper [25] at the mean field level the moments c n show slightly 
more structure. The rise in C2 is somewhat sharper, the peak in C4 is about 5 % higher. 
In summary, however, the corrections induced so far by corrections involving (<£* — $) 
around the mean fields are small. Pionic fluctuations, to be discussed in Sec. 5, tend to be 
more important. In presently available lattice results [9], these latter effects are however 
suppressed by the relatively large pion masses. 

The ratio of the moments C4 and C2 has been discussed [7] as a suitable indicator of 
fluctuations close to the phase transition. The quantity of interest here is the cumulant 
ratio i?4 2 defined in [7] and given as R\ 2 = I2C4/C2. The PNJL model calculation for 
this ratio is shown in Fig. 8. The dashed curve is found in the mean field limit with 
($*) = (<J>) which suppresses one of the two Polyakov loop degrees of freedom. The solid 
curve is computed with inclusion of corrections beyond mean field and demonstrates the 
role of the non-zero (<£* — <&). At temperatures below T c one reaches R\ 2 = 9, the value 
characteristic of a hadronic resonance gas [7]. 
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Figure 8: The cumulant ratio R\ 2 from 
the PNJL model in and beyond mean field 
approximation in comparison with lattice 
QCD results, with C2 and C4 as given in 
Ref. [9]. 



Figure 9: Contributions to the conformal 
measure from quarks and Polyakov loop, as 
well as the total PNJL interaction measure 
for N t = 2. The N { = 2 + 1 lattice QCD re- 
sults [10] (with iV r = 8 for p4-improved and 
asqtad action) are shown for orientation. 



4.4 Interaction measure 

The PNJL results for the interaction or conformal measure e — 3p are illustrated in Fig. 9. 
The total interaction measure normalised to T 4 is split into quark and Polyakov loop parts. 
Note the sensitive balance between quark quasiparticle and Polyakov loop contributions 
to e — 3p close to T c . In pure gauge QCD (or with infinitely heavy quarks) the Polyakov 
loop interaction measure is positive throughout. The presence of light quarks and their 
dynamical coupling to the Polyakov loop changes their pattern significantly. The Polyakov 
loop parts of the presure itself that determines the dashed curve in Fig. 9, is found to 
be consistent with calculations reported in Ref. [27]. For orientation, the total PNJL 
interaction measure (with Nf = 2) is shown in Fig. 9 along with recent TVf = 2 + 1 lattice 
QCD results [10]. 



4.5 Speed of sound 

In Fig. 10 the squared speed of sound in units of the speed of light is plotted as a solid 
line. The speed of sound v s is defined by 



dp 




on 








~d~e 


s 


df 


vl 


dT 2 


V 



(4.9) 



where the denominator is the specific heat capacity cy. The dashed line in Fig. 10 gives 
the size of the ratio of pressure and energy density, |. In the panel to the left where 
the quantities are plotted at vanishing chemical potential fi = both graphs show a 
pronounced dip near the crossover transition temperature. In the panel to the right the 
same situation is plotted at a quark chemical potential close to the chemical potential of 
the critical point fx < ^ C rit.- 
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Figure 10: The speed of sound (solid) and the ratio of pressure over energy density (dashed) 
at vanishing chemical potential as function of temperature (left panel). The right panel 
shows the same quantities at a quark chemical potential slightly less than the one at the 
critical point (fx = 0.3 GeV < fi crit ~ 0.31 GeV). 



5 Dynamical fluctuations 

So far the formalism presented has been focused on the treatment of fluctuations around 
the mean fields, averaged over space and (Euclidean) time. The homogeneous, constant 
Polyakov loop field and its corrections beyond mean field fall in this category. 

In this section we consider mesonic excitations and their propagation (i. e. -^--correc- 
tions). The NJL framework is well suited to incorporate such effects. The NJL model 
features a dynamical mechanism which produces spontaneous chiral symmetry breaking 
and, at the same time, generates the pion as a Goldstone boson in the pseudoscalar quark- 
antiquark channel, together with a massive scalar (sigma) boson. The thermodynamics 
of these modes and their changing spectral properties have been subject of several NJL 
model calculations in the past [32-34]. 

With increasing temperature, the mass of the pion is still protected by its Goldstone 
boson nature, whereas the sigma mass drops until at T ~ T c it becomes degenerate with 
the pion, signalling restoration of chiral symmetry in its Wigner-Weyl realisation. For 
T > T c , the 7r and a masses jointly increase quite rapidly while at the same time their 
widths for decay into qq grow continuously. This implies that at temperatures exceeding 
T c both 7r and a modes become thermodynamically irrelevant while correlated quark- 
antiquark pairs carrying the quantum numbers of ir and a can still be active above T c . 
One therefore expects that the corrections to the pressure from propagating pions and 
sigmas should be concentrated around T c . These mesonic modes are colour singlets 9 . 
Thus their statistical weight is much smaller than the weight of the deconfined quark 
quasiparticles. 

9 Colour octet quark-antiquark modes turn out to be heavy and far removed from the spectrum of active 
degrees of freedom. 
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5.1 Meson propagators in the PNJL model 

We start from the derivation of mesonic propagators in the PNJL model as performed, for 
example, in [35]. We calculate the momentum dependent propagator 



k 



d 2 s } 



bos 



(5.1) 



where £ = 9 — 9mf now stands for the pion field or for the deviation of the sigma field 
from its expectation value. Note that the functional trace in the formula for tS D os ensures 
momentum conservation, such that the sum of the momentum arguments in the denom- 
inator always vanishes. The calculation can be done numerically as it was done in the 
previous section. Alternatively, we use an analytic approach as follows. Recall some useful 
formulae also exploited in Refs. [6,9]: 



8 In det M 



dx 



= tr 



AT 



_! dM 
dx 



and 



dM' 1 „„ ,dM % , , 

— " ~ ~dx • (5 ' 2 > 



with M an invertible matrix and is the component-wise derivative of this matrix. 
Applying this to the PNJL action iS D os in (2.8) and neglecting the potential terms for the 
moment we find 



as, 



bos 



V 



86 



d 3 p 
(2^ 



Tr 



S {iu n ,p ; 9) — 



(5.3) 



where S -1 (iu n ,p; 9) denotes the inverse quark propagator with emphasis on the fact that 
the quark propagates in the mesonic background field 9. This formalism makes it possible 
to calculate derivatives with respect to bosonic fields (say 9k) that have not been explicitly 
included in the action, as long as it is ensured that the model does not produce finite 
vacuum expectation values for these particular fields. Not having a vacuum expectation 
value is equivalent to the fact, that the mean field equations corresponding to these fields 
are satisfied for a vanishing field, i. e. that 



as, 



bos 



89 



. 



(5.4) 



e k =o 



All we need to know is the constant matrix 9 g 1 . This matrix involves the Dirac, colour 
and flavour structure of a quark-antiquark pair (or a quark-quark pair) that couples to 
the bosonic field 9k, i. e. it is determined by the quantum numbers of 9k- The mean field 
equation is fulfilled if the trace in Eq. (5.3) vanishes for the given Dirac, colour and flavour 
structure. For the pion field this is true as long as there is no pion condensate. 10 The 
condensate corresponding to the sigma, namely the chiral condensate, figures explicitly in 
the action and is therefore included in the quark propagator. 

In the case of the pion and sigma propagators the functional derivative in (5.1) produces 
exactly the trace over Dirac, colour and flavour structures known from RPA calculations 
[32-35]. We adopt the definition of the quark distribution functions f£ and and the 
separation of the emerging integral into the contributions I\ and I2 as given in Ref. [35]. 



°The mean field equation is satisfied as the flavour-trace trf [ 1 n] = with i — 1,2,3 vanishes. 
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In the treatment of the thermodynamics we have modified the cutoff prescription of the 
standard NJL model, such that non-divergent integrals are integrated over the whole 
quark-momentum range, while only divergent integrals are regularised by the usual NJL 
three-momentum cutoff. The separation of finite and divergent contributions is defined 
such that the model reproduces the classical limit at high temperatures, i. e. the Stefan- 
Boltzmann limit. As a downside, for consistency all newly appearing integrals have to be 
treated in the same manner, which leads to slightly different results from those given in 
Ref. [35]. 



5.2 Mesonic corrections to the pressure 

Once the meson propagators are given, it is possible to evaluate the contribution to the 
pressure from mesons propagating in the heat bath using RPA methods. Applying Bethe- 
Salpeter (RPA) equations generates spectral functions 

/ - ^ GImn M (w, q; T) 

Pm(uj, q; T) = 



(l-GRen M ) 2 + (GImn M ) 2 
with the thermal quark-antiquark polarisation function 

U M (uJ,q; T) = J J^3 Tv ^MS(iu n + fj,,p)TuS(i(u n -u)+fj,,p-q) 

where the sum is taken over the Matsubara frequencies u n = (2n + 1)itT. Here Tm is a 
Dirac, flavour and colour representation of a meson current labelled M. In this work we 
only focus on the pseudoscalar isovector channel (i. e. pionic excitations) and the scalar 
isoscalar channel. S(iu n ,p) = — ^ ™Jz+ m z denotes the quark quasiparticle propagator 
with f = iuj n j - 7 ■ p. 

The pressure below T c is essentially generated by the pion pole with its almost tempera- 
ture independent position. Therefore the calculated pressure below T c basically represents 
the one of a pion gas with fixed (temperature independent) mass. Fig. 11 shows, as exam- 
ples, the spectral functions for the pion and sigma modes at threshold temperature T t h r 
where the breakup into a quark-antiquark pair occurs. This threshold temperature is at 
about 1.1 T c . At this point the ir and a spectral functions are still distinguishable (left 
panel of Fig. 11), whereas they coincide (right panel) at temperatures well above threshold 
where ir-a degeneracy indicates restoration of chiral symmetry in its Wigner-Weyl realisa- 
tion. Their width is a measure of the decay of the (increasingly massive) pionic and sigma 
modes into (light) deconfined quark-antiquark pairs at temperatures above T c . 

The resonant interaction of instable mesons with the quark sea above T c produces 
an additional pressure contribution. This contribution is not part of the quark pressure 
previously calculated in Hartree-Fock approximation. The meson decay products form 
rings of RPA chains. Such kind of pressure contributions are investigated in Ref. [36] and 
calculated performing the ring sum. However, below T c the NJL model does not handle the 
mesonic degrees of freedom properly. In the hadronic phase the coupling of mesonic modes 
to the quark-antiquark continuum is suppressed by confinement, whereas pm receives 
contributions from decays into qq even below T c . This unphysical feature persists [35] 
in the PNJL generalisation of the NJL approach. Moreover, the non-renormalisability 
of the NJL model requires to introduce further subtractions when following the lines of 
Ref. [36] . To avoid such arbitrariness and unphysical features we ignore the decay of meson 
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Figure 11: The spectral functions pm = (i^GRcUM^+fhimUM^ taken at q = for pion 
and sigma at T « T t h r (left) and at T > T t hr (right). 



modes into (/(/-pairs altogether when calculating an estimate for the meson contributions 
to the pressure: 

m=v S Tln(1 - e ~ Eq/T) + B[T) ' (5 - 5) 

where v is the statistical weight of the corresponding meson species, E q = q 2 + wip ole (T) 

with m po i e (T) the temperature dependent pion and sigma pole mass determined by 1 — 
GReLT = 0. Furthermore B(T) is an appropriately chosen vacuum energy constant en- 
suring thermodynamic consistency. B(T) is fixed such that the temperature dependence 
of the pole mass m po i e (T) is compensated on differentiating 0, with respect to the temper- 
ature T. This implies that the inclusion of B(T) ensures that dQ,/dm po \ e \ T = 0. 

In Fig. 12 the calculated pressure of 7r 0,± and sigma modes are compared with the 
quark Hartree-Fock pressure and the result for the overall pressure of Hartree-Fock plus 
RPA is plotted. For comparison the pressure of a Bose gas with three internal degrees of 
freedom is indicated by the thin solid line. Below the crossover temperature T c one can 
clearly identify the pion gas contribution resulting from the RPA calculations. Once the 
meson masses reach the scale of the NJL cutoff A the used approximation breaks down. 
The inversion of the scale hierarchy appears at temperatures of about 1.3 T c . 

For larger current quark masses the meson gas contributions and correlations are re- 
duced. This effect is illustrated by Fig. 13 where the pressure of the PNJL model is plotted 
using an increased current quark mass leading to an unphysically heavy pion. Thus for 
heavy pions the agreement with lattice data observed in a previous publication [25] re- 
mains. This agreement is also confirmed by calculations in a non-local PNJL framework 
[37] which does not suffer from cutoff artefacts. 

The mesonic contribution to the interaction measure is rather modest. The interaction 
measure already shown in Fig. 9 is replotted in Fig. 14 including mesonic contributions. 

6 Conclusions and outlook 

The PNJL model as an approximation to QCD thermodynamics picks up on two basic 
properties of low-energy QCD: spontaneous chiral symmetry breaking and confinement. 
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Figure 12: The pressure contribution origi- 
nating from pion modes, sigma modes and 
from quarks in Hartree-Fock approximation 
(dotted). The thin solid line represents the 
pressure of a gas of bosons with three inter- 
nal degrees of freedom and a constant mass 
m = m^iT = 0). 



Figure 13: Same as Fig. 12, but with higher 
current quark mass tuq = 50 MeV =>■ m n = 
421 MeV (compared to m = 5.5 MeV =¥■ 
ma = 139 MeV in Fig. 12). The pressure 
of the boson gas (thin solid line) was now 
plotted using the heavier pion mass. 
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In this work some of the existing calculations [19, 25, 26] have been extended in several 
directions. We have reviewed the expectation values of the Polyakov loop and its complex 
conjugate, the phase diagram, the moments of the pressure and the speed of sound in a 
framework beyond mean field theory. While the phase diagram does not show significant 
changes when improving the mean field approximation, the moments of the pressure and 
the speed of sound show quantitative differences on the order of 5%. In general the 
structures observed become more articulate. In the case of the Polyakov loop and its 
complex conjugate the corrections cause qualitative differences. While the Polyakov loop 
and its complex conjugate are equal at mean field level in the present approach, the 
corrections beyond mean field generate the split of the two expectation values ($) and ($>*) 
at non-zero quark chemical potential. The numerical results show that the corrections are 
largest in the vicinity of phase transitions or rapid crossovers. This comes as no surprise 
as it is the transitional region between two phases where we expect large fluctuations. 

The degrees of freedom that govern the low temperature regime, primarily the pions, 
produce significant corrections to the pressure only in the regime below the critical tem- 
perature T c where constituent quarks are frozen and confined. As soon as the pressure of 
quark degrees of freedom starts to rise at the chiral and deconfinement crossover, mesonic 
pressure contributions become comparatively small. 

The good agreement of the PNJL model at mean field level with lattice calculations 
remains in the presence of the mesonic corrections calculated in this work. The most 
prominent feature of the pressure, namely the steep rise near the critical temperature, is 
only slightly modified by the corrections due to dynamic fluctuations of pions and sigma 
mesons. Below the quark-antiquark threshold the pressure generated by pion fluctuations 
is basically the pressure of a free pion gas. As low temperatures are difficult to access 
by lattice QCD, the pressure in typical lattice calculations is usually normalised to zero 
at some finite temperature below T c . This might explain why the pressure of the PNJL 
model including mesonic corrections is slightly higher than the pressure resulting from 
lattice calculations [9]. Due to these normalisation issues the comparison suffers from 
this uncertainty, A (p/T 4 ) = j>/T 4 „_„ — p/T 4 L, n , which in turn depends on the 

I 1 — 1 norm ' ^ *U 

normalisation temperature T norm and additionally on the realized pion mass. For large 
pion masses (see Fig. 13) this correction is small maintaining the good agreement between 
PNJL and lattice results. For small pion masses the pressure contribution from pion 
modes is almost flat in the temperature region T ~ m^, such that the pressures from 
lattice and PNJL calculations mainly differ by a shift in p/T 4 . Shifting the lattice data to 
higher values of p/T 4 indeed reduces the difference between Stefan-Boltzmann limit and 
lattice data for the pressure at high temperatures around 2-3 T c and above, improving 
the agreement between lattice results and PNJL. Even when taking into account these 
issues in the comparison of PNJL and lattice results, we conclude that there exists a good 
qualitative and quantitative agreement of these two approaches. 
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A Detailed derivation of corrections to mean fields 



A.l Expansion of the effective action 

This appendix displays some technical details concerning the treatment of fluctuation 
corrections beyond mean field approximation in the PNJL model (cf. Sec. 3). 

In the following we denote by 6 = (Oi) the set of fields (a, A, fa, fa) which operate as 
bosonic degrees of freedom in the effective action 5b os of Eqs. (2.8) and (2.11). Further- 
more, let 6*o = ((c) j (A) , (fa) Q j (fa)o) be the set of mean field (expectation) values of 
these quantities, and introduce deviations from the mean fields by £ = = 9 — 9q. 

A frequently used procedure that we follow here, is to expand the effective action in 
powers of £ around a properly chosen mean field configuration. The Gaussian part of such 
an expansion of the path integral can be handled analytically. In Sec. 3 the mean field 
approximation has been defined such that the (formally) complex action 5b s produces, to 
this leading order, a real- valued thermodynamical potential (or pressure), Omf = R< e [^o]> 
subject to the mean field equations (2.16). The expansion of 5t>os is then of the generic 
form 

5bos = ^(^MF + ^ (1) -e+^^ (2) ■£•••) , (A.l) 

where we have introduced the notations a ■ b = ^ a, b\ and a ■ A ■ b = Ylij a i A%i 1 
with summations extending over all bosonic degrees of freedom. The expansion (A.l) is 
performed such that the path integral is optimally approximated. This is achieved when 
the perturbative terms in the expansion of the action are maximally suppressed. With 
the thermodynamic weight e~ s £ C this approximation is optimal near the maximum of 
|e [. The equations to determine #0 are the mean field equations (2.16) (also used in 
[25])- 

Given the expansion (A.l) in terms of the £ fields, thermal expectation values incor- 
porate fluctuations around the mean field configuration 6mf = $o- We refer to these 
corrections as "fluctuations" even if the fields themselves (such as the Polyakov loop field 
variables fa and fa) are constant in space and time. 

A perturbative approach is now used to calculate corrections to the mean field solu- 
tions. The action iSbos is split into "large" and "small" parts, 5bos = So + <Si, as follows: 
the "large" part So incorporates the leading mean field terms plus the additional Gaussian 
part of 0(£ 2 ) in Eq. (A.l): 

So = ^(Re[fto] + ^ (2) -£) , (A.2) 

while S\ deals with the remaining pieces, in particular with the non- vanishing Im[fio]- 
The leading correction of this sort is the term 5Si = ^ru/ 1 ) • £. In the present context we 
truncate Eq. (A.l) as it stands and keep only this term in Si, for the moment. 
The thermal expectation values of a given quantity /(£) is proportional to 

J Vi f(0 e" 5 "" = jvi f(0 e- 5 « e-* , (A.3) 

where, for fields constant in space-time, the path integral reduces to 

J d£ f(0 e- s ^ e- s ^ = J d£ f(0 e^>«) e~ lk < , (A.4) 



22 



= -ImcjW . (A.5) 
(i. e. about 6 = #mf) in powers of £ involves 

J de r e - s »(« e -*< = (id k r z (k)\ k= v luiw(1) , (a. 6 ) 

where we have introduced the generating function Zo(k) = J d£ e~ s °^ e~ ik '^. Each 
power of idk evidently produces a factor y. At the same time, performing this derivative 
explicitly on Zo(k), with <So(£) specified in Eq. (A. 2), produces a factor 

(A.7) 

which is independent of y. 

Hence there are two small quantities at hand to establish a perturbative expansion: y 
and 5. The smallness of y is given here as we are interested in the thermodynamic limit. 
The size of 5, however, is controlled by the action itself. Whether the expansion in 5 is 
justified or not depends on the model and must be examined accordingly. The explicit 
calculations presented in the main body of this work shows that in the present version of 
the PNJL model the expansion in 5 is indeed a good approximation. 

We are now in a position to write down the thermal expectation value of a generic 
function / as an expansion in powers of y and S. We proceed here with establishing 
Feynman diagrams for this perturbative approach. We write generically 

Z = jfJ Vte~ s ^ = Ifjvt fji (-«Si)< e- 5 «. (A.8) 

If corrections to the partition function of the PNJL model are to be calculated, the So 
part of the action only comprises zeroth and second order terms, while the "small" part Si 
is identified with all other orders. The first order term acts as a source term. We establish 
the following Feynman rules: 



with 

k = 

iT 

A perturbative expansion of /(£) about £ = 
integrals of the form 



S = l v 



(2) 



— 1 



■ k 



,(2) 



dS hos 
d 3 S hos 



+ 



I 
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j 






k 



d 2 Sbos 

d 4 £bos 



(A.9) 



In perturbation theory it can be shown that only connected diagrams contribute to 
the partition function, i. e. 

oo (- oo ^ 

Zi = (e- 5l )o = E n <(-^i) Z >o = exp £ ~ { ((-^Oc , (A.10) 

1=0 ' ln=l ' J 
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where (• • -) denotes the expectation value with respect to the unperturbed action, and 
(• • -) 0c is the expectation value of the connected diagrams with respect to this unperturbed 
action. Note that here the corrections depicted by the Feynman diagrams are corrections 
to the negative action, —S, as the partition function was defined by Z = e~ Sctt -. The 
corrections therefore need to be subtracted from the mean field result of the action 5mf- 
For the thermal expectation values of / we write 

oo ^ 

(f) = (fe' S % = z ^j^(f(-S 1 ) l ) . (A.ll) 
Here each term under the sum can be written in terms of connected expectation values 

(f = E ai !a 2 !(2!)^...'(a n !)(n!^m! ^~ S ^ ^^ol ' ' ' 

a\,a,2--- ,a n ,m=0 

■■■(/(-SHoAi, (A.12) 
where v = a± + 2a>2 + • • • + na n + m. Substituting back in Eq. (A.ll) gives 

C 00 ^ 00 

(/ T*>o = e XP E ((-*r) 0c x Yl ~ x (f (-^i) m ) 0c • (A-13) 

(n=l ' ) m=0 

Using Eq. (A. 10) we find the final result 

00 1 

(/) = (/ e- 5l )o = E ~ { (/ (" 5 i) n )oc • (A-14) 

n=0 

In terms of Feynman diagrams Eq. (A. 14) can be translated into all those connected 
diagrams that contain exactly one insertion coming from the function /. The Feynman 
rules for the insertions of / are 



df_ 



d 3 f 



d2 f 

o = - 

j k dijd^k 

I 



m 



d 4 f 



i atjd&db 



k 



dijd^kdiidir, 



(A.15) 



What is needed to use these rules systematically is a scheme that orders all possible 
diagrams according to their importance in powers of the small parameters y and 5. The 
lowest order corrections in y and 5 are shown in Table 1. 

A useful consistency check is to verify that the thermal expectation values are now 
closer to the properties of an order parameter than the mean field result. In other words: 
we examine whether the thermodynamic potential Q is a Landau effective action minimised 
with respect to (a) , (A) , (<J>) , (<£*) using Eq. (A. 14) for the expectation values. The 
analysis below is done for the lowest order terms, a = and (3 = 0, 1. We start from the 
form also used for the numerical calculations, presented below Eq. (3.1), and differentiate 
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/(#mf) 
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a = 1 






2 1 |x * > X 










a = 2 


*e 
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4 v — r x / j \ — 

i>KX>- 





Table 1: The Feynman graphs contributing to (/), ordered in i^y) and 5^ with multi- 
plicity factors. 
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= 1 


a = 






a = 1 


\ x 2 o — o 


I X 4 o — o — x + |x2 ^ — x 



Table 2: The lowest order Feynman graphs contributing to xjj- The vertices depicted 
as a circle are now the contributions of g, defined analogously to the contributions of / 
in Eq. (A. 15). The prefactors are the product of the multiplicity factors of the original 
Feynman graphs and the factors arising from the differentiation. 
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with respect to the expectation values (9) = ((a) , (A) , (<fr) , (<£*)) T . To orders a = and 
(3 = 0, 1 we find that (0) = 9q + 50, where (50 is given by 
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(A.16) 



(which is Eq. (3.2) with /(6>) = #j ). After some calculation we arrive at the lowest order 
term in (3 
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8n 



8 2 n 
~~bW 



8n 

~~89 



+ higher orders, (A.17) 



which is of order (3 = 2, i. e. the self consistency equations are satisfied to the order we have 
been working in. As a consequence the corrections necessary to account for the fermion 
sign problem do not modify the mean field equations. A backward reaction on the mean 
field equations does not occur at this level of the approximation, which focuses on the 
evaluation of corrections to the effective (Polyakov loop) potential. Additional effects of 
pionic and scalar quark- antiquark modes, as considered in Sec. 5, do in principle have 
backward effects on the mean field equation. 

The formalism allows to determine susceptibilities involving a quantity g, \g = [V((9 2 }~ 
(g) 2 )] 1 / 2 . All that needs to be done is to apply the previously developed formalism 
to the function g 2 . In Table 1 the Feynman rules and multiplicity factors are writ- 
ten down for the evaluation of (/). In a second step / is replaced by g 2 . In this step 
the product rule of differentiation has to be applied producing additional prefactors. In 
this procedure it will happen that vertices of / with m = 2, 3, . . . or more legs will 
split into two vertices with mi + rri2 = m legs. The lowest orders of the expression, 
Eq. (4.1), are shown in Table 2. The contributions of order (|J) cancel. In this frame- 
work susceptibilities scale with V? as expected. Additionally, it becomes obvious from 
Table 2 that there are no mean field contributions to susceptibilities in the sense that 

0.9 ~ (5')mf) 2 )mf = (# 2 )mf ~ (#)mf = 9mf ~ 9mf = °- In the framework of mean field 
calculations, susceptibilities are usually evaluated by inverting the second derivative of the 
mean field action with respect to the fields. This is seen in the present framework as well: 
the entry for a = 1 and (3 = in Table 2 produces exactly this expression. 
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